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Abstract 

We show how Jarzynski relation can be exploited to analyze the nature of order-disorder and 
a bifurcation type dynamical transition in terms of a response function derived on the basis of 
work distribution over non-equilibrium paths between two thermalized states. The validity of 
the response function extends over linear as well as nonlinear regime and far from equilibrium 
situations. 
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I. INTRODUCTION 



Advancement of micromanipulation techniques has opened up the new possibility of esti- 
mating equilibrium thermodynamic quantities from non-equilibrium measurements in recent 
years^H^. The experiments with pulling forces of piconewton magnitude on RNA molecules^ 
and dragging colloidal particles in a fluicP^ have made it possible to measure the probability 
distribution of work exerted on a system. Since the relaxation time of the system concerned 
in these experiments is too long compared to the time scale over which the non-equilibrium 
measurements are made, work fluctuations over non-equilibrium irreversible paths play an 
important role in estimating the free energy difference between the thermalized initial and 
final states of the system. This is reflected in celebrated Jarzynski relational which relates 
the work distribution and the Helmholtz free energy difference AF as 

e -PAF = ( e -W) ( L1 ) 

where the averaging on the exponential function of the work variable W has been carried 
out with work distribution P(W), (3 being 1/ksT with ks and T denoting the Boltzmann 
constant and temperature, respectively. Central to these studies are several results, e. g., 
the so-called fluctuation theorems expressing transient violation of second law of thermody- 
namics where the systems in question are driven arbitrarily far from equilibrium. 

The probability distribution of work P(W, t) done on a system by manipulating an exter- 
nal agency or force is followed over irreversible paths for a given protocol. As free energy is 
the key quantity for carrying thermodynamic information and for description of transition 
between different states, it is apparent that, by virtue of Jarzynski equality^ it is possible 
to probe the signature of dynamical transition in the behaviour of work distribution. The 
problem is non-trivial since for a nonequilibrium system, in general, free energy functional 
can not be defined. This, in consequence, raises difficulty in constructing a response function 
directly in terms of free energy and its derivatives with respect to suitable parameters of the 
system in the spirit of what is done in equilibrium phase transition. This difficulty, however, 
can be overcome with the help of Jarzynski equality by probing work fluctuations on the 
system even far from equilibrium. The object of the present paper is to examine this issue 
and to look for the response function which is characteristic of the dynamical system under 
study and is derivable from work distribution. In what follows we use a model potential to 
explore the order-disorder transition as well as a transition arising out of bifurcation of the 
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response function in terms of the time evolution of work distribution function. The dynam- 
ics is described by a Langevin equation which appropriately takes care of the variation of 
the potential parameter in inducing transition from an unimodal to bimodal character* 13 * 14 !. 
The implication of the results are elaborated. 



II. THE MODEL AND THE EQUILIBRIUM DESCRIPTION 

Consider a Brownian particle in a thermal bath at temperature T and subjected to an 
external potential force. The governing Langevin equation is given by 

x = p (2.1a) 
p = - lp -V'{x,t)+T(t) (2.1b) 

where V(x,t) is the external potential and 7 is the dissipation constant, x and p denote 
the coordinate and the momentum of the particle, respectively. Thermal fluctuations T(t) 
of the bath are modeled by Gaussian, zero mean and delta correlated noise 

(r(t)) = (2.2a) 
(r(t)r(t')) = 2D6{t-t') (2.2b) 

where D = ^ksT is the strength of the noise. We make use of potential V(x, t) of the 
following form 

V{x,t) = \(t)a }x 2 + h)x A (2.3) 

where a, a , b are the potential parameters and X(t) is a switching parameter. We first 
allow the system to reach an equilibrium with the heat bath at temperature T and then 
switch on the parameter, infinitely slowly from an initial state A = to a final state A = 1. 
The model has been explored in the literature under diverse conditions 13 > 14 . At the initial 
state the potential is bistable with two minima at x = ±^a/b and one maximum at x — 0. 
The corresponding distribution function is bimodal in position coordinate (x) 



P {x,p) ~ £(p)exp 



l/2ax 2 - 1/4 bx 



(2.4) 



D 

where £(p) = exp (— p 2 /2KBT)At the final state, A = 1, for a > ao the potential is also 
bistable with two shifted minima at x — ±^ (a — a )/b and one maximum at x — 0. For 



a = ao, the potential contains only x 4 term and the corresponding distribution function is 
unimodal as follows; 

"-1/4 bx 4 ' 

P/(z,p)~£(p)exp -±- (2.5) 

So the switching process induces a transition of the distribution function from bimodal(A = 
0) to unimodal one(A = 1). This is depicted in Fig.l(a,b). Fig.l(a) also includes the case 
for a < ao for A = 1. The potential well is steeper compared to previous case (V(x) = bx 4 ). 
The corresponding probability distribution is also shown in Fig. 1(b). 
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FIG. 1: A schematic illustration for switching of (a) the potential energy for different values of the 
parameter A around the transition point for (a — A(i)ao) = a corresponding to Vq (A = 0, solid 
line), (a — A(t)ao) = a/2 corresponding to Vt (A = 1/2,, dashed line, (a — X(t)ao) = corresponding 
to Vf (A = 1, dotted line), and dash-dot line present the case (a — A(t)ao) < (b) distribution 
function for the corresponding cases. In all the illustrations a = a®. 

As both the initial and final states are in equilibrium the free energies of these states are 
given by 

F = -k B T In Q ; F f = -k B T In Q f 

respectively, where Q denotes the partition function. These are given by the following 
expressions 

Qo 



/+oo 
dp dx exp 
- 



Qf 



oo 
oo 



(1/2) p 2 - (1/2) ax 2 + (1/4) bx 4 



knT 



dp dx exp 



(1/2) p 2 - (1/2) (a - a )x 2 + (1/4) bx 
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(2.6) 
(2.7) 
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respectively. By numerical integration one can easily find out the change in free energy for 
arbitrary values of a . 



AF = F f — Fa = —knT In 



Or 



(2.8) 



The linear coefficient ao serves here as the constant parameter of the system, ao governs the 
nature of "phase" as a region of space in which the free energy function is analytical and 
continuous. The dynamical transition is associated with the crossing of boundary between 
the two regions. To understand the nature of transition it is therefore worthwhile to look 
for any discontinuity or irregularity of the first and second derivatives with respect to the 
control parameter at the boundary. With this in mind we, in Fig 2(a), present the variation 
of free energy change as a function of ao at different temperature. In the spirit of traditional 
way of deriving thermodynamic information and characterizing equilibrium phase transition, 
one may now define an extensive variable, an order parameter type quantity as, — (^§^ j 
analogous to magnetization M associated with the free energy of a magnetic system as 
M = — (f§) T (therefore role of magnetic field H is played by a ). Pushing the analogy a 
bit further we define a response function as a second order derivative with respect to a as 

d 2 AF , x 

* = -in < 2 - 9 > 

Clearly the above expression is analogous to the magnetic susceptibility which measures 
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FIG. 2: (a) Variation of free energy change with ao at different temperature, (b)variation of 
with ao at different temperature, for the parameter set a = 0.5, b = 0.005, Kb = 1. 



the variation of magnetization due to the change of the external field. The magnetic sus- 
ceptibility is related to the second order derivative of free energy as \ — ~ {w^- ^ n Fig-2 
(b) we depict the first order derivative of free energy change with respect to ao at different 
temperature. While no discontinuity or singularity is observed in the curve, the behaviour 




FIG. 3: Response function vs ao plot at different temperature for the parameter set a = 0.5, b = 
0.005, K B = 1 (using Eq.(2.9)). 

of response function (2.9) as shown in Fig. 3 is markedly different. The response function 
gradually increases to a maximum followed by a sharp fall at around a = a. This type of 
variation of response function with ao illustrates two types of dynamical transition^ 2 , (i) 
The first one corresponds to a transition from an ordered state to a dis-ordered state for 
AV = ksT. The barrier height of the potential (AV) is (a — ao) 2 /Ab. So with increasing 
value of a (ao > 0) the barrier height gradually decreases. At ao = a — 2y/bkBT, the barrier 
height is double to average kinetic energy of the particle. The response function^) at this 
point shows a maximum, presumably indicating the existence of a disordering transition. 
The existence of the transition point can be physically explained as follows: Whenever the 
value of average kinetic energy crosses half of the barrier height the particles start moving 
more randomly between two wells in a way as if they (relatively large number of the parti- 
cles) feel no barrier between them. This is due to the fact that in this temperature domain 
the distribution of kinetic energy gains a broad range. Here lies a valid reason for existence 
a relation between barrier height and kinetic energy in the transition region. As expected 
with the increase in temperature the transition point shifts towards the origin, (ii) The 
second type of transition is due to a change of distribution function as reflected in Fig. 3 
by a very sharp change (bifurcation) in the response function at ao = a. This transition is 
intrinsically different in nature from the previous one which can be controlled by adjusting 
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the temperature of the system. In the later case, however, the temperature has no significant 
influence. 



III. JARZYNSKI RELATION AND WORK DISTRIBUTION FOR TRANSITION 
IN NONEQUILIBRIUM SYSTEM 



For infinitely slow switching(A = to A = 1) of the potential parameter(Eq.(2.3)), 
the system remains in quasistatic equilibrium with the reservoir throughout the switching 
process and the total work performed on the system will be equal the Helmholtz free energy 
difference between initial and final state states (W = AF). If the switching process occurs 
with a finite rate, the total work spent in changing the state of a system will depend on 
the microscopic initial conditions of the system and the reservoir and obey the following 
inequality 

(W) > AF 



In this case, the time evolution of the system of interest as described earlier(Eq.(2.1a,b)) is 
governed by a stochastic phase space trajectory, which depends on the externally imposed 
time dependence of the switching parameter X(t). For the present purpose we consider a 
constant switching rate, A = l/t/(t/is the final time starting from t = 0). With this time 
dependence of switching parameter the total work performed along one particular trajectory 
(z(t)) up to time t is given by 

W[z{t),t] = £dt (^)f^(O) (3-1) 
The equation for time evolution for work is given by 

Although the equation of motion for W does not have an explicit dependence on noise it is 
stochastic through its dependence on phase space variables. Before we proceed to analyze 
the essential features of dynamical transition in terms of the work done we first numerically 
simulate the stochastic coupled equations Eq( 2.1a| b) along with the equation for work (3.2) 



simultaneously using standard Heun's algorithm and calculate the distribution of work, av- 
erage work, higher moments of work. This allows us to calculate the free energy change 
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FIG. 4: Plot of work distribution function for different values of switching rate for the parameter 
set a = 0.5, b = 0.005, T = 1.0, Kb = 1 and 7 = 2.5. 

from Jarzynski relation as a non-equilibrium estimate. We use in our numerical simulation 



very small time step(At) of 0.01 for numerical integration has been used. For the initial 
conditions we have assumed that at t — all the particles are in the potential minimum at 
x = \[ajb with zero velocity. In our simulation, we first allow the system to equilibrate with 
the reservoir to smooth out the effects due to the influence of initial conditions and transient 
processes. After the equilibration process we switch on the parameter (A(t)). In calculation 
of average and higher moments the averaging is done over 20,000 trajectories. The paramet- 
ric dependence of stochastic trajectory takes care of energy balance-^ for Langevin dynamics 
when work is calculated in terms of Eq.(3.1). 

We first calculate the work distribution, average work and work fluctuation for different 
switching rate. This is shown in Fig.4 and Table-I. As revealed by Fig. 4 the distribution 
is nearly Gaussian for low switching rate centering around (W), which is larger than free 
energy change ( AF) (The vertical dotted line does not indicate the center of the distri- 
bution but shown only to highlight the portion at which W matches AF). With decreas- 
ing switching rate the center of the distribution is shifted towards the line of free energy 
change, AF = 11.862536 (average work done tends to be equal to the free energy change 
along with the decrease of width of distribution). The distribution function P(W) tends to 
P(W) 5((W) — AF) . In Fig. 5 we have depicted the work distribution function along 



a slowly varying time dependent quantity X(t), with switching rate A ~ 10 5 — 10 6 . A 
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FIG. 5: In this figure we have depicted work distribution function for different switching rate with 
corresponding Gaussian fitting curve. The parameter set: a = 0.5, b = 0.005, T = 1.0, Kb = 
1 and 7 = 2.5. 

TABLE I: The table presents nonequilibrium estimates of average work and width of work distri- 
bution for different values of switching rate for a = 0.5, b = 0.005, T = 1.0, Kb = 1 and 7 = 2.5. 



Switching rate 


(W) 


AW 2 

(= (W 2 ) - (W) 2 ) 


9.375 x 10~ 5 


11.889384 


0.049791 


3.125 x 10~ 5 


11.878023 


0.032513 


9.375 x 10~ 6 


11.870123 


0.017320 


3.125 x 10~ 6 


11.863661 


0.002563 



with corresponding Gaussian fitting curve for relatively larger switching rate (A > 10~ 5 ). 
From Fig. 5 it is apparent that the work distribution function is not a Gaussian for arbi- 
trary switching rate. One observes that the higher order cumulants are important. (The 
distribution will be a Gaussian only if the switching speed is very slow for which the higher 
order cumulants effectively zero). We present Table-I depicting the variation of average work 
and width of the work distribution with switching rate. To verify the Jarzynski relation 
we numerically estimate the quantity (e~ w ^ kBT ). For a comparison of this result with the 
analytically calculated free energy we present a data Table-II for different values of a . The 
calculated free energy change using Jarzynski relation 

AF = -k B T In ((e~ w/kBT )) (3.3) 
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TABLE II: The table presents a comparison of free energy change calculated analytically (equilib- 
rium estimate) and numerically using Jarzynski relation (non-equilibrium estimate) for different 
values of ao- The other parameter set for this data table are a = 0.5, b = 0.005, T = 1.5, Kb = 
1, 7 = 2.5, and Switching rate = 5 x 10 -6 . 



a 


AF 

(using Eq.(2.8), 
equilibrium 

estimate ) 


(W) 


AW 2 


AF (using Jarzynski 
relation Eq.(3.3 ) 
non-equilibrium 
estimate) 


(fluctuation-dissipation 
estimate) 


0.46 


11.452822 


11.456438 


0.011571 


11.452579 


11.452581 


0.47 


11.522062 


11.525216 


0.0119362 


11.521237 


11.521238 


0.48 


11.58799 


11.592910 


0.012375 


11.588780 


11.588785 


0.49 


11.650821 


11.655346 


0.012639 


11.651129 


11.651133 


0.50 


11.710751 


11.714050 


0.012955 


11.709732 


11.707931 


0.51 


11.767964 


11.772241 


0.013103 


11.767872 


11.767873 


0.52 


11.822633 


11.827798 


0.013666 


11.823242 


11.823243 


0.53 


11.874918 


11.878662 


0.014022 


11.873984 


11.873988 


0.54 


11.924969 


11.929216 


0.014020 


11.924543 


11.924542 



matches almost exactly (absolute difference is less then 0.01 %). As revealed by the data 
set of Table-II the following relation for the free energy change with work holds very good 
as a fluctuation-dissipation estimate 

AF = (W) - {(W 2 ) - (WO 2 } (3.4) 

the dissipative work being equal to the width of the distribution of work, 

w diss = ^ {(W 2 ) - (W) 2 } (3.5) 

A final remark on Table-II may be in order. Since the work fluctuations have been determined 
with the help of Langevin dynamics with a Gaussian noise, fluctuation-dissipation estimate 
in Eq.(3.4) matches well with free energy change. For a different protocol for following the 
trajectory or for higher switching rate the deviation from linear response is expected where 
non-Gaussian distributions of work make their presence felt. The latter aspect is evident in 
Fig.5. 
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IV. STUDY OF DYNAMICAL TRANSITION BY WORK FLUCTUATION AND 
RESPONSE FUNCTION 



We recall that a system out of equilibrium can not, in principle, be described by free 
energy and therefore no free energy functional or partition function methodology can be 
applied for non-equilibrium system for classification of dynamical transition. Guided by a 
close analogy with equilibrium phase transition in Sec. II we have identified an appropriate 
response function \1/ for description of dynamical transition. Analysis of Sec. Ill on the other 
hand suggests that by virtue of Jarzynski relation one can relate the free energy change 
(rather than free energy itself) between two thermodynamic states to microscopic work fluc- 
tuation for non-equilibrium paths connecting these states. The question is, can we bypass 
the description based on free energy change as done in Sec. II to compute the response func- 
tion \l/ directly from work fluctuation and recover the features of non-equilibrium dynamical 
transition. To this end we now proceed to calculate the response function \1/ in terms of the 



work performed using Eq.(3.4) and Eq.(3.5) as follows; 



, r _ d 2 AF _ d 2 (W) | d 2 w dlss 
dal dal da^ 



From Eq.(3.1) and (2.3) the work performed for the switching process along a particular 
trajectory is given by 

W[z(t)) = jf ' dt (I) \a x 2 (4.2) 

As revealed by the above expression the work performed is a linear function of ao only if the 
process x(t) is independent of control parameter (ao). In this situation we have 

d(W) = (WO _ d 2 (W) = Q 
da a ' dal 



With the help of the Eq.(4.3) the response function can be expressed in a simplified form as 



1 (W 2 ) - (W ) 2 
a 2 2k B T 



The response function is thus related to dissipative work. The above expression for response 
function has a close similarity to other linear response functions, as for example, heat capac- 
ity C v = ^ E 2kgT^ ' ma g ne tic susceptibility x — ^2A^t^ > where E and M denote internal 
energy and magnetization, respectively. A difference between linear response functions and 
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FIG. 6: The solid and dotted lines present the variation of response function without dissipative 
work as a function of ao for different values of temperature. The lines containing hollow and solid 
circles present the variation of response function with inclusion of dissipative work as a function of 
ao for different values of temperature. The parameter set used a = 0.5, b = 0.005, 7 = 2.5 Kb = 1 
. The inset plot presents the variation of work as a function ao for the same parameter set as in 
the main figure but for T = 1.0, Switching rate = 10~ 5 . 

however, is noteworthy 11 . Since E and M are well-defined equilibrium properties of a 
macroscopic system in the thermodynamic limit while W by its very nature corresponds a 
quantity denned for non-equilibrium paths, heat capacity and susceptibility are typical static 
response function in contrast to the non-equilibrium response In the present problem 
the work, however, is not a linear function of ao- This is due to the fact that a bears an 
implicit dependence on x(t). The non-equilibrium estimation of (W) shows that, in general, 
it is a nonlinear function of a , as depicted in the inset of plot Fig. 6. (W) behaves linearly 
only for small values of a . The second derivative of the average work also predicts two 
transition points as discussed in Sec. II. But the transition points are slightly shifted from 
original position of the transition point (this is depicted in Fig. 6). It should however be 
emphasized that although the second order variation of average work predicts only the ap- 
proximate transition point, the inclusion of the contribution due to dissipative work gives 
an accurate estimate of the transition points as shown in Fig. 6 (the plots containing solid 
and hollow circle). 

In order to check the dependence of response function on switching rate we present the 
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FIG. 7: The sub-figures (a, c, e) present the variation of response function \t with ao and the 
subfigures (b, d, f) present the variation of $ = — 9 with ao depicting dynamical transition for 
different switching rate. The parameter set used a = 0.5, b = 0.005, 7 = 2.5 Kb = 1. 

variation of response function as a function of ao for different switching rate in the Fig. 7 
(a,c,e). One observes that fluctuations become larger for higher switching rate. It may, 
however, be checked that the nature of the variation of response function with ao remains 
same. Although the nature of variation of the response function indicates two types of 
dynamical transition for arbitrary switching rate, but it is difficult to point out the transition 
points for higher switching rate due to fluctuating nature of ^. Moreover in order to figure 
out the relative contribution of second derivative of average work, i. e., — i^r 1 and that of 
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dissipative work m res P onse function we present the sub figures of Fig. 7 as a pair 

(a,b), (c,d) and (e,f). While Figs.7(a,c,e) represent variation of response function including 
dissipative work, Figs.7(b,d,f) represent the cases without dissipative work. From these 
figures it is clear that one can identify the dynamical transitions by calculating only second 
derivative of average work, i. e., 9 !£P ■ Thus the nature of the variation of response functions 
derived from work fluctuations are practically independent of any pulling speed which is a 
correct reflection of Fig. 3. The response function(4.1) or (4.4) is completely determined by 
the stochastic equations (2.1) and (3.1) of the dynamical system independent of free energy 
description of the system in the thermodynamic limit. 



V. CONCLUSION 



Analysis of transition between two equilibrium states is traditionally based on free energy 
change in a system with respect to a relevant parameter and linear response function is 
related to fluctuation of the order parameter around equilibrium. In view of the fact that 
free energy remains undefined for non-equilibrium systems, any analysis of response function 
on the basis of similar argument is untenable. However as Jarzynski equality is related to 
free energy change with work fluctuations during the passage of the system over many 
non-equilibrium paths when the parameter of the potential is varied, it is possible to look 
for the signature of dynamical transition in the behaviour of work fluctuations. The key 
quantity is the response function, which exhibits a characteristic behaviour of the system 
itself with no limitation imposed by near-equilibrium condition. Based on a model system 
we have examined two types of transition, e., g., order- disorder type and a bifurcation type 
which can be differentiated by their thermal behaviour around the transition points. As the 
Jarzynski relation and the related fluctuation theorems are valid even far from equilibrium 
situations, the response functions derived from work fluctuations, we believe, have a wider 
range of applicability, i., e., beyond linear regime and the results obtained for this simple 
system may be extended to explore more complex issues. 
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